Three-dimensional MHD simulation of expanding magnetic flux ropes 
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Three-dimensional, time-dependent numerical simulations of the dynamics of magnetic flux ropes 
are presented. The simulations are targeted towards an experiment previously conducted at CalTech 
(Bellan, P. M. and J. F. Hansen, Phys. Plasmas, 5, 1991 (1998)) which aimed at simulating Solar 
prominence eruptions in the laboratory. The plasma dynamics is described by ideal MHD using 
different models for the evolution of the mass density. Key features of the reported experimental 
observations like pinching of the current loop, its expansion and distortion into helical shape are 
reproduced in the numerical simulations. Details of the final structure depend on the choice of a 
specific model for the mass density. 

PACS numbers: 52.30.-q, 52.65.-y, 52.30.Cv 



I. INTRODUCTION 

Magnetic fields play an important role for the struc- 
ture and the dynamical behavior of the solar atmo- 
sphere. Well-known examples for structural features 
found on a variety of spatial scales are helmet streamers 
sunspots, coronal loops and filaments, while flares, 
loop eruptions and coronal mass ejections are dynamical 
phenomena related to magnetic energy release ||. The 
importance of the magnetic field structure and its evo- 
lution has prompted a lot of theoretical investigation in 
the past. For example, Torok et al. || use numerical 
simulations to study a scenario for loop eruptions due to 
kink mode instabilities, using the loop model by Titov 
and Demoulin 

An entirely different approach to study the evolution 
and creation of magnetic signatures found in the solar 
corona is the use of laboratory experiments. Bellan and 
Hansen report from an experiment which is considered a 
laboratory model for Solar prominence eruptions [H, Q , 
and a similar experiment ( "FlareLab" ) has recently been 
set up at Ruhr-Universitat Bochum. In the latter, the 
long term goal is to study a variety of magnetic field con- 
figurations and to employ extended plasma diagnostics. 

Fig. [T] shows a crude sketch of these two experiments: 
A horse shoe magnet is mounted below the bottom plate 
of a vacuum chamber, producing an arc-shaped field in- 
side the vessel. Two electrodes in the vicinity of the 
magnet poles are connected to a capacitor bank. In 
the experiment, hydrogen gas is puffed into the cham- 
ber shortly before the capacitor voltage is connected to 
the electrodes. A plasma arc is created by ionisation 
of the hydrogen, causing the voltage to break down and 
the capacitors to be discharged. The evolution of this 
arc is then followed by means of a fast camera in order 
to study the dynamics that follow from the internal mag- 
netic forces. Images from different stages of the discharge 
process as published in [5] indicate that the arc pinches 
in cross-section and expands as a whole on the scale of 
a few microseconds. At a later stage, it gets further de- 
formed and assumes a helical shape. A corresponding 
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FIG. 1: Sketch of the initial experimental setup used by Bel- 
lan [g| and in "FlareLab" at Bochum Universitat. 



photograph recently taken from "FlareLab" is shown in 

Fig. m 

Later, this basic experimental setup has been modified 
by Bellan et al. to include an additional magnetic field 
component with field lines that cover the entire filament 
structure at a right angle. Using this modified setup, it 
was demonstrated that the "strapping" field component 
can inhibit the arc expansion Q. Also, the interaction 
between two current arcs has been studied in an other 
extension of the basic configuration [f| . 

This experimental approach offers an attractive way 
to gain insight into the dynamical behavior of magnetic 
structures, because it opens a way to analyze details of 
the dynamics through the measurement of physical quan- 
tities like the magnetic field, which are not directly ac- 
cessible in the solar context. Apart from this, the pos- 
sibility to control and selectively modify parameters like 
boundary- and initial conditions allows systematic inves- 
tigations and the test of theoretical models or comparison 
with computer simulations. However, a number of ques- 
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II. MODEL AND NUMERICAL TREATMENT 



FIG. 2: Image taken from the "FlareLab" experiment at 7ps 
after ignition. The plasma arc has already expanded to ap- 
proximately three times its initial radius and helical distor- 
tions are clearly visible. 



tions arise concerning the interpretation of the experi- 
mental observations and its relevance for solar physics: 
First of all, it is not clear if the model of magneto- 
hydrodynamics (MHD), which seems appropriate in the 
solar context and has been adopted for the interpreta- 
tion of experimental findings [|[, can be applied to the 
experiment without limitations. While the plasma arc 
itself may be highly ionized, there is no indicaction avail- 
able from the experiment concerning the electron density 
off the regions of highest luminosity. Therefore, it is un- 
clear if, for instance, the "frozen-flux" principle applies 
in those surrounding regions and what the consequences 
of any "non-ideal" behavior would be. 



We employ the model of ideal MHD, 
-V ■ (pv) + S 



dp 
dt 
dv 
dt 
dB 
~dt 
V x B 



(1) 



-(tf.VW+^^-^ + j/Aff (2) 
P P 



V x [v x B 
Poj 



(3) 
(4) 



Here, p, v, B, p and j denote the mass density, plasma 
bulk velocity, magnetic field, thermal pressure and elec- 
tric current density, respectively. A homogeneous kine- 
matic viscosity, v, is included for reasons of numerical 
stability. In |T]), a yet unspecified source term S allows 
the model to deviate from plasma mass conservation. As 
described below in detail, this term is used in order to 
implement different models for the plasma density in our 
simulations, namely i) mass conservation, S — 0, ii) a 
homogeneous density, p = po iii) a fixed Alfven velocity, 
p oc \B\ 2 and iv) a crude model for ionization and recom- 
bination (case 4 below). These models are only partly 
motivated on physical grounds, but we note that Eq. (|2|) 
is written in non-conservative form and hence remains 
consistent with |T]) even for S ^ if it is assumed that 
plasma created according to S has the same velocity as 
the existing plasma population. For an ionization term, 
this is the case if the neutral gas component moves at 
the same velocity v as the charged species. To close Eqs. 
dXI) — d3J) 3 we assume an isothermal plasma, p = Tp, in case 
i) (i.e. mass conservation) and p = in all other cases. 

All quantities are normalized to a typical value of the 
magnetic field strength, Bq, a length scale, Lo, and the 
Alfven velocity, va = Bo/ ^popo with a typical value for 
the mass density, po. A normalization time scale follows 
as t = Lq/va, and the relationship of these scales with 
the experimental setup are described below. 



Initial conditions 



In an attempt to obtain a better understanding for the 
dynamical evolution of the plasma arc, we have carried 
out a number of computations which aim at simulating 
the plasma dynamics in the experiment. Using the model 
of ideal MHD, we can address the question whether ex- 
perimental results in fact are compatible with that model, 
and to which extent it is capable of reproducing the ob- 
served structures. In addition, detailed data is available 
about relevant quantities like the magnetic field, elec- 
tric current and plasma flow. It should be stressed that 
these simulations are directed at the laboratory exper- 
iment alone, and no attempt is made at this point to 
draw conclusions for prominence dynamics in the solar 
atmosphere. 



At t = 0, the plasma is assumed to be at rest, v = 0, 
with a prescribed density distribution according to one 
of the specific density models. Hence, we make no at- 
tempt to describe the ionization stage of the experiment, 
but rather start with a configuration that contains highly 
conducting plasma in the entire domain. 

A cartesian coordinate system is used in which the 
plasma chamber wall that contains the two electrodes 
lies in the plane z — 0, and the electrodes themselves are 
located at ±Re y . 

The initial magnetic field consists of two parts: The 
contribution from the horse shoe magnet in the experi- 
ment is modeled by the field B^x) created by two virtual 
magnetic dipoles located outside the domain at positions 
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x± = ±Re y — Re z , i.e. "below" the electrodes. They 
are assumed to carry dipole momenta m± = ±me z with 
m chosen such that the field at the electrode positions is 
normalized, \Bd(±Re y )\ = 1. It follows that 

^ _ 3 (mi ■ (x - Xi) (x - Xi)) - rhi\x - Xi\ 2 



A second initial magnetic field component, B c (x), 
mimics the field related to the plasma current in the 
early stage of the discharge. A quantitative model for 
B c is constructed by means of a vector potential A c [x) 
with B c (x) = Vx A c (x), and the choice of A c (x) is moti- 
vated by the assumptions that i) the initial current is ap- 
proximately localized around the half-circle in the plane 
x = which connects the electrodes, and ii) the direc- 
tion of A c coincides with the direction of B^ because the 
ring current will roughly follow those field lines. With 

ed{x) :— given from the dipole field model above, 

we choose 

A c (x) = A e -(l-l-«) 2 /* 2 e - 2 / s2 e d {x) ( 6 ) 

The exponentials localize the magnetic structure on a 
scale of 8 = 0.625 around the ring of radius R in the y-z- 
plane. The amplitude Aq is determined from the condi- 
tion that the maximum of \B C \ in the bottom plane z = 
is max 2= o(|i? c |) = 3, which roughly corresponds to the 
ratio of the current-induced magnetic field to the horse- 
shoe magnetic field magnitude taken at the position of 
the electrodes as estimated for the experimental setup 
(see below). 

At this point, we would like to stress that the resulting 
initial magnetic field, B(x, t = 0) = 5d(f) + Vx A c (x), is 
not force-free. Even with a (small) thermal pressure term 
added, there is no force equilibrium at the start of the 
simulation, which is in accordance with the experiment. 
The aim here is, similar to the experiment, to investigate 
dynamical properties of the resulting evolution. 



B. Density models 

In the next section, we will describe results obtained 
from four simulation runs, all of which use the same ini- 
tial magnetic field as described above, but differ in the 
treatment of the mass transport, i.e. in the specific real- 
ization of Eq. (fl}: 

Case 1 - Mass conservation: Here, Eq. (fT]) is used 
as a continuity equation with S = 0, describing mass 
conservation. More specifically, we use a homogeneous 
mass density p — 1 as initial condition and specify the 
temperature as T = 1. This model can be seen as the 
simplest choice possible that accounts for a consistent 
mass transport during the dynamical evolution of the 



system. Realizing that \B\ « 3 in the current filament, 
the resulting local plasma-/? is 2/\B\ 2 w 0.2 <C 1. 

Case 2 - Fixed density: We keep p = 1 fixed 
throughout the simulations, i.e. Eq. (fT]) is abandoned. 
As a consequence, sound waves are eliminated from the 
dynamics. This model is used in order to get an indica- 
tion of the influence of the mass transport and pressure 
included in the previous case on the evolution. 

Case 3 - Fixed AlfVen velocity: In this case, 
the mass density is continuously adjusted such that the 
Alfven velocity is constant throughout the domain, i.e. 
p(x,t) oc \B(x, t)\ 2 . The interest for this study stems 
from the fact that the evolution might be treated as 
quasi-static |5j, which means that the Alfven crossing 
times are small compared to the global evolution time 
scale. Setting va to a constant value results in a homo- 
geneous communication of Alfvenic disturbances. 

Case 4 - Ionization/recombination model: Here, 
we implemented a simple model for the ionization by the 
electric current and chose the source term as 

S = T i j 2 + T r (p-p ). (7) 

The coefficients are set to I\ = 0.5, T r — 5 and po = 1 so 
that the time scales of ionization and recombination are 
comparable to the Alfvenic and convection time scales. 



C. Numerical implementation 

Equations (fT])-©, with a specific density model, are 
discrctized by finite-differences on a cartesian grid and 
integrated as an initial value problem using a third 
order Runge-Kutta scheme. The numerical box cov- 
ers [-20,20] x [-20,20] x [0,40] in the x-, y- and 
z— direction, respectively. Boundary conditions at z = 0, 
i.e. the "bottom" plane, are such that v = and B± 
is linearly extrapolated. We use outflow boundary con- 
ditions, i.e. the ghost cell values are set to the first cell 
values within the domain. However, the simulated cur- 
rent filaments are well separated from those boundaries 
so that these conditions have no significant influence on 
the results. 

In order to obtain sufficient spatial resolution of the 
current filament dynamics, we carried out mesh-adaptive 
computations with local refinement using the simulation 
code "racoon" [9(. Grid blocks of 16 3 cells each are recur- 
sively refined up to a local resolution equivalent to 1024 3 , 
using a refinement criterion that compares the electric 
current density \j\ with a critical value j cr u that is de- 
termined from the existing local resolution. A typical 
grid layout is depicted in Fig. [3] 
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FIG. 3: The full computational domain with the block distri- 
bution and the line-of-sight integrated electric current density, 
viewed parallel to the :r-axis. 



D. Scaling of experimental parameters 

The model described above involves a number of nor- 
malization values for the MHD quantities. In order to 
allow a direct comparison with the experiments of Bellan 
et al. and "FlareLab," we relate these values to physical 
parameters given in [B[ and estimates from the "Flare- 
Lab" experiment (not yet published). 

With Lq — 16 mm, the electrode spacing is 5Lo = 8cm 
which corresponds to the experimental value, and the en- 
tire domain covered in the computation is 64 x 64 x 64cm 
which roughly corresponds to the dimensions of the 
plasma vessel. Taking Bq — 0.3 T, which matches the 
horse shoe magnetic field at the electrodes in the "Flare- 
Lab" experiment, and assuming a hydrogen plasma of 
no w 10 20 m~ 3 , results in an Alfven velocity of Va = 
6-10 5 m/s and an Alfven crossing time of to ss 2.6 -10~ 8 s. 
This value is small compared to the macroscopic evolu- 
tion time scale of rs 1/xs reported for the experiment, but 
two things should be noted here: First, from experiment 
the detailed density distribution can only be estimated as 
the experiments still lack suitable diagnostics to measure 
this value. More importantly, the magnetic field strength 
drops off drastically with increasing distance from the 
magnet poles/current filament, and the filament length 
increases in time from its initial value of f=a 10 Lq. Hence, 
the true Alfvcnic travel time along the filament from one 
electrode to the other will be longer than to and might 
reach the overall evolution time scale. 

As for the magnitude of the initial magnetic field cre- 



ated by the plasma current, B c , we assume the entire dis- 
charge current of 7 « 50 kA to be located in a channel of 
radius 5 — 0.625Lo = 10mm. By Ampere's law, this will 
create an azimuthal magnetic field of B v = ^,qI/2tt5 — 
1 T « 3.3 Bo on the channel surface. Hence, we chose 
the magnitude of A c in Eq. ([6|) such that |£? c | = 3 close 
to the electrode. 

The viscosity in ((2]) is included for numerical stability 
of the simulations. It was necessary to use a normalized 
value of v = 0.01 in cases 1 and 2 and an even larger 
value of v = 0.05 in cases 3 and 4. These values might 
correspond to a significantly higher viscosity than the ac- 
tual physical one of the experiment. However, not only is 
the latter largely unknown by value, but it is also ques- 
tionable if the viscous term in the momentum equation 
provides a realistic description at all. Hence, we treat 
the viscosity as a purely numerical parameter necessary 
for stabilizing the computations. 



III. SIMULATION RESULTS 

The observations in the experiment are made with an 
optical CCD camera. It is obviously difficult, if not im- 
possible, to unambiguously relate the luminosities of the 
images to plasma quantities, but we see it as a reasonably 
assumption that there exists a close correlation between 
the electric current density and the luminosity. There- 
fore, the simulation results presented here focus mainly 
on the current density evolution, and volume renderings 
of \j\ are shown in the plots discussed below for compari- 
son with the images taken by the laboratory camera. Fur- 
ther, a view perpendicular to the current arc is chosen, 
again to fit the camera's perspective in the "FlareLab" 
set-up. 



Fig. [4] shows the spatial distributions of \j\ and p at 
three different times of simulation case 1, i.e. the run im- 
plementing mass conservation. Essential dynamical fea- 
tures can be observed on different time scales: Within 
the first Alfven time, the current arc pinches to approx- 
imately half of its initial diameter. This effect can be 
determined from the plots taken at t = 5.1: The current 
is more localized than in the initial configuration, and 
the mass density is enhanced in the arc center as a con- 
sequence of the compressional plasma motion. Accord- 
ingly, the density is reduced around the arc with mini- 
mum values of ~ 0.53. This deviation of mass density 
distribution from the initial homogeneous value corre- 
sponds to the formation of pressure gradients due to the 
assumption of isothermal plasma with normalized tem- 
perature T — 1. It should be noted that the pinching is 
less pronounced close to the electrodes where the dipole 
magnetic field component is comparable in magnitude to 
the non-potential B c . Apart from this pinching, the arc 
as a whole has undergone an expansion from its original 
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FIG. 4: Electric current density \j\ (left) and mass density p (right) of simulation run 1 at times t = 5.1 (top), t = 12.3 (center) 
and t = 43.2 (bottom). 
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major radius R = 2.5: The apex has risen toz~ 4.5 and 
the curvature in the upper part of the arc is reduced com- 
pared to the original half-circle shape. Also, close to the 
electrodes, a more pronounced deformation of the current 
channel into helical shape is observed. Both, the expan- 
sion and the foot point deformation occur on the scale of 
several Alfven times, i.e. slightly slower than the pinch- 
ing. In the following evolution, the current structure and 
the mass density perturbation continue to expand (resp. 
times t = 12.2 and 43.3), causing the loop apex to rise 
to about z = 10, i.e. four times its initial height. In the 
course of this expansion, the tendency of helical defor- 
mation manifests itself on the entire structure and grows 
in amplitude. Eventually, the upper part of the loop is 
entirely rotated out of the plane x = and intersects that 
plane at right angle with its apex (time 43.2). Later, the 
expansion slows down and finally stalls shortly after the 
last frame shown in Fig. [U After that, the structure 
actually starts to shrink slightly but finally comes to rest 
in a configuration close to the last one shown. 

The qualitative interpretation of this sequence in terms 
of ideal MHD appears to be straight-forward: Recalling 
that the initial configuration is not in force balance, the 
current will pinch towards a force balance perpendicular 
to the center of the current tube. With the original diam- 
eter of 26 = 1.25 and an amplitude of the current-induced 
field of \B C \ — 3 and density p ~ 1, this process occurs on 
a time scale of « 0.1. The expansion of the arc can be in- 
terpreted as the response to the well-known "hoop force" 
caused by the arc curvature. As the arc follows this force, 
the curvature and the corresponding current density are 
reduced and the expansion slows down. In the vicinity of 
the foot points, the line- tying condition and the zero ve- 
locity boundary condition prevent the arc from following 
the expansion in the horizontal directions, and hence the 
arc looses its circular shape with its upper part merely 
rising upwards. Finally, we interpret the formation of he- 
lices at the foot points as kink modes and the fact that 
these modes develop fastest at the electrodes as a conse- 
quence of the different local Alvfen velocities: While the 
magnetic field is dominated by B c on the entire arc, the 
mass density is considerably reduced to values around 
0.5 close to the electrodes in the early stage dynamics 
(cf. p at t = 5.1 in Fig. Hence, the Alfvenic kink dy- 
namics will be twice as fast as in the apex, where p « 2 
as a consequence of the early pinching. This estimate is 
consistent with the fact that kinks become apparent at 
t ~ 10 at the apex, while they form already around t ss 5 
at the electrodes. At later stages, the entire structure re- 
laxes into an approximately force-free state in which the 
internal twist of the magnetic field has been converted 
into a large-scale writhe. 

Using the results of case 1 from above as a reference, 
we describe in the following the significant differences 
observed in the remaining simulation runs. 



Fig. [5] displays the current structure of case 2 ("fixed 
density") at different times. While the overall evolution 
is similar to case 1, some consequences of the constant 
mass density, and hence the absence of pressure forces, 
are readily visible. First of all, the maximum current 
density during the pinching phase is significantly larger 
than in case 1 (« 45 compared to « 10 at time t « 1.5), 
apparently because of the absence of restoring thermal 
pressure forces in case 2. Another consequence of the ho- 
mogeneous mass density is that the kink formation in the 
apex region occurs already at t w 5, at the same time as 
it is observed close to the foot points. This homogeneity 
along the arc slightly alters the overall structure forma- 
tion: While case 1 displayed an evolution from the foot 
points to the apex and from short wavelengths to larger 
scales, the development in case 2 is almost uniform along 
the arc. Small scale kinks occur early and grow in am- 
plitude until the perturbation of the arc becomes almost 
stationary in shape, and only the slow expansion of the 
structure can be seen in the simulation. In particular, 
the pronounced rotation of the structure as a whole as 
seen in case 1 is absent here. Rather, the late stage can 
still be interpreted as a large current arc lying in the 
x = 0-plane and continuing to expand. 



Fig. [5] shows that the evolution of the current struc- 
ture in case 3 ( "fixed Alfven velocity" ) resembles closely 
the results of case 2. The upper part of the arc flat- 
tens and forms kinks, however, the formation of kinks is 
constrained to the upper part of the arc and the "legs" 
remain almost unperturbed. 



Finally, the inclusion of a simple recombina- 
tion/ionization model in case 4 again leads to a 
slight differences in the resulting current system (Fig. 
[7J, but the overall shape is still comparable to the 
previous cases, in particular to case 2. 

Comparison with Experiment and Discussion 

As noticed before, the experimental observations and 
the numerical simulations can only be compared on a 
qualitative level, not least because of the missing quan- 
titative data from the experiment: It is not clear in 
which way the photographic luminosity originating from 
the neutral hydrogen emissions is related to the physical 
plasma parameters and fields. Based on the assumption 
that the electric current is the dominant source for ion- 
ization, and thus the plasma density and temperature 
will be well correlated to the current density, we choose 
to compare the laboratory images to \j\ obtained from 
the numerical simulation. 

Further fundamental differences between the experi- 
mental regime and the simulation model must be kept in 
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FIG. 5: Electric current density \j\ of simulation run 2 at times t = 1.5, t = 5.1, t — 18.7 and t = 49.6. 



mind: First of all, the model of ideal MHD is an ideal- 
ized approach to the physics of the plasma discharge, in 
particular with respect to the degree of ionization, the 
ionization process during ignition, the role of collisions 
and possible kinetic or at least multi-fluid effects. In ad- 
dition, the treatment of the boundaries and the coupling 
of the electrodes to the plasma is highly idealized. An 
other difference between model and experiment is that 
the total discharge current increases in time during the 
arc evolution, while the simulation boundary conditions 
imply a constant, prescribed current through the lower 



plane. Concerning the time scales, we note that a di- 
rect identification of simulation time with observed time 
scales is difficult because the fundamental plasma param- 
eters of the experiment are currently only estimated. 

Keeping these limitations in mind, the computations 
reproduce key features of the laboratory results like the 
arc expansion and the kink formation. For example, 
by comparing Fig. [2] with the distributions of \j\ from 
the four simulation runs, the qualitative agreement be- 
comes evident. The moderate modifications of the den- 
sity model that we employed in the four different runs 
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FIG. 6: Electric current density \j\ of simulation run 3 at times t = 2.2, 7.2 and 9.6, respectively. 



lead to slightly different details in the arc dynamics as 
manifested, e.g. in the location, number and intensity of 
kinks. Asking which of these models can best reproduce 
the experiment, we observe that a characteristic feature 
there is the fact that the central region of the arc con- 
sistently gets bent downwards (see Figj2]). This fact has 
been reported by Bellan et al. @ as well. In the simu- 
lations, we observe the same behavior in run 3 (constant 
Alfven velocity) and we explain it by the fact that the 
mass density oc B 2 assumes its largest values in the apex 
of the plasma arc. The reason for this is that the initial 



dipole field Bd is weakest there which leads to a relatively 
strong pinching effect and hence large values of B. The 
large mass density close to the apex causes that part of 
the arc to lag behind during the arc expansion, leading 
to the characteristic dip. 

We repeated the simulation of case 3 with an increased 
value for the viscosity [y = 0.1 instead of the previous 
v = 0.05) and achieved an even better agreement with 
the experimental pictures (Fig. [8] a) . With this modifi- 
cation, the overall structure of the arc appears smoother 
compared to the previous run 3. The data of this high vis- 
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FIG. 7: Electric current density \j\ (left) and mass density (right) of simulation run 4, time t = 19.9. 



cosity run has been used to produce a perspective view 
of the arc structure, showing magnetic field lines and 
electric current lines (Fig. [8]b). In this view, the over- 
all structure of the arc as a helically distorted current 
channel becomes much more evident than from the axes- 
parallel views and shows an excellent agreement also with 
an image published in Q . 

Taking these findings together, we conclude that the 
present simulations represent a successful initial step to 
model the overall dynamics of the experiment in the 
framework of ideal MHD. Despite the fact that the 
present models contain a number of ad-hoc assumptions 
(like a specific density model, the presence of viscosity 
and a simple isothermal energy transport equation), the 
results are remarkably encouraging. For the future, we 
expect improved experimental diagnostics to give impor- 



tant additional information that in turn might lead to 
a more adequate theoretical and numerical description 
and hence to an even better agreement between both ap- 
proaches. 
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FIG. 8: a) \j\ from run similar to case 3 but with v = 0.1 at time t — 18.0. b) Perspective view of selected magnetic field lines 
(green), electric current lines (red) and \j\ (grey shade) from the same run. 



